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Over the past few decades there has been an enormous 
growth in the numbers of papers devoted to experimental 
and theoretical aspects of anomalous diffusion [H, H L The 
landmark review by Metzler and Klafter in 2000 l| has 
been particularly influential, promoting the description 
of anomalous diffusion within the framework of continu- 
ous time random walks (CTRWs) and fractional calcu- 
lus. There are now numerous applications utilizing this 
approach in physics, chemistry, biology and finance [ij. 

A central theoretical result in this research was the 
derivation 0, 0] of a fractional Fokker-Planck (Smolu- 
chowksi) equation Q 



D 

dt "0^* 



P{x,t) (1) 



to describe the evolution of the probability density 
function P{x, t) for subdiffusion in an external space- 
dependent force field F{x). In this equation. 
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is the Ricmann-Liouvillc fractional derivative, IS a 
fractional diffusion coefficient, rja = (/3Ka)~^ is a frac- 
tional friction coefficient, and /3 is the inverse tem- 
perature /csT. The fractional Fokker-Planck equation, 
Eq.([T]), was derived from the continuous time random 
walk model of Montroll and Weiss 0, with power law 
waiting times 0, 0] ■ 

More recently a modified fractional Fokker-Planck 
equation. 
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'P{x,t), (3) 



was derived from power law waiting time CTRWs, us- 
ing a generalized master equation [7[, for subdiffusion in 



a time-dependent force field F{t). The modified frac- 
tional Fokker-Planck equation, Eq.Q, was also derived 
for subdiffusion in dichotomously alternating force fields 
Bill, F{x)i{t) with ^(t) ±1, but a fractional Fokker- 
Planck equation to model subdiffusion in general space- 
and-time-dependent force fields F(x, t) has remained elu- 
sive [S- 11|. On the other hand, a subordinated stochastic 
Langevin equation has been proposed for modelling sub- 
diffusion in space-and-time-dependent force fields [loj . 
More recently [l2[, in the case of time-dependent forces, 
the moments of the stochastic process defined by this 
stochastic Langevin equation were shown to coincide with 
the moments of the modified fractional Fokker-Planck 
equation, Eq.®. 

There have been numerous papers on fractional 
Fokker-Planck equations in recent years relying on ad 



3nt years relyi: 



hoc or phenomenological models 

In this letter we derive the fractional Fokker-Planck 
equation. 
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from power law waiting time CTRWs, using a generalized 
master equation, for subdiffusion in a space-and-time- 
dependent force field F{x,t). This fractional Fokker- 
Planck equation is shown to be formally equivalent to 
the subordinated stochastic Langevin equation in [loj for 
space-and-time-dependent forces. We also show that the 
original fractional Fokker-Planck equation, Eq.(IT]), gen- 
eralized by replacing F(x) with F(x, t), can be recovered 
from power law waiting time CTRWs in an ad-hoc gen- 
eralization of the CTRW particle balance equation if the 
diffusing particles respond to the force field at the start 
of the waiting time prior to jumping. These derivations, 
and further extensions to include reactions, are described 
in greater detail, for chemotactic forcing, in a related 
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publication [15| . 

Our starting point is the generalized master equation 
approach developed in 0, This approach utilizes 

two balance conditions. The balance equation for the 
concentration of particles, ni{t), at the site i and time t 



dnjjt) 
dt 



(5) 



where jf are the gain (+) and loss (-) fluxes at the site 
i. The second balance equation is a conservation equa- 
tion for the arriving flux of particles at the point i. In 
general, to allow for biased CTRWs in a spacc-and-time- 
dcpendcnt force field, we write 

J+it) = prix,^ut)Jr_,{t) + pi{x,+,,t)J-^,{t), (6) 

where Pr{x, t) andp;(a;, t) arc the probabilities of jumping 
from X to the adjacent grid point, to the right and left 
directions respectively. The two balance equations can 
be combined to yield 



dni(t) 
dt 



■Pr{x^-i,t)J^_^{t) +pi{xi+i,t)J.^^^{t) - J, [t). 

(7) 

For CTRWs with a waiting time probability density func- 
tion il}{t) the loss flux at site i is from those particles that 
were originally at i at t = and wait until time t to leave, 
and those particles that arrived at an earlier time t and 
wait until time t to leave, hence 11611 



J-{t) = ^(<)n,(0) + i>{t-t )J+{t ) dt . (8) 



We can combine Eq. (O and Eq. ([5]) to obtain 



jr{t) ^^P{t)n,{0) + J 4'{t~t) 
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dt , 
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'(10) 

where the hat denotes a Laplace transform with respect 
to time and s is the Laplace transform variable. This 
simplifies further to 



and then 

jr(s) = i){s)n,{0) + i^{s) J' (s) + sn,{s) - n^{0) 



(11) 



where $(s) is the Laplace transform of the survival prob- 
ability 



$(t) = j ^{t')dt'. 



(12) 



In the CTRW model, subdiffusion originates from a 
heavy-tailed waiting-time density with long-time be- 
haviour [l[ 



T \ T 



(13) 



where a is the anomalous exponent, r is the characteristic 
waiting-time, and k is a dimensionless constant. Using a 
Tauberian (Abelian) theorem for small s 17 



Aa - 



(14) 



where Aa 



/ta(l — a) 



In the special case of a Mittag- 



Leffier waiting time density 18| the ratio, Eq. (fTl)) . is ex- 
act and Aa = 1. The loss fiux can now be obtained 
by using the ratio, Ea (jl4[) . in Eg. pip and inverting the 
Laplace transform. Noting that the Laplace Transform 
of a Ricmann-Liouvillc fractional derivative of order a, 
where < a < 1, is given by [l9| 



d'^fjt) 
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this yields 



Aa di-"n,(i) 



dt 
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(15) 
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where we have assumed that the last term in Eq. <\15\i is 
zero. Using this result in Eq.Q yields. 



dnt{t) _ Aa 



dt 



— <^Prix.,-i,t)- 
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+Pl{Xi+i,t)- 



dt^-°^ 
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dt 
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The jump probabilities are biased by the external 
space-and-timc-depcndent force. Here we consider (near 
thermodynamic equilibrium) Boltzmann weights with 



and 



PriXi,t) = C 



Pl{Xi,t) = C 



exp(-/3F(x»+i,t)) 
eM-PV{x,,t)) ' 

exp(-/3y(a:,-i,t)) 
exp(-/3y(x,,t)) 



(18) 



(19) 



The jump probabilities are determined at the end of the 
waiting time, when the particle must jump, so that 



Pr{Xi,t) +pi{Xi,t) = 1, 

which defines C, and then 



Pl{Xi,t) -pr{Xi,t) 



-l3V{xi-i,t) 



-fiV(xi+i,t) 



-pv{xi^ut) _|_ f,-PV(xi+i,t) ■ 



(20) 



(21) 
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The spatial continuum limit of Ea. p7|) can now be 
obtained by setting Xi = x and Xi±i = x ± Ax and 
carrying out Taylor scries expansions in x. Retaining 
terms to order (Aa;)^ and using Eq. ipiH) yields 



dn{x, t) 
dt 



^ |ax^ [{pi{x,t) -pr{x,t))^Dl-"n{x,t 
Ax^ 



+- 



-,D^-°'n{x,t) 



2 dx^° * 

The Taylor series expansion of Eq. ((21]) yields 



(22) 



Piix, t) - Prix, t) « /3Ax^4^ + OiAx'), (23) 



and then 

dn{x, t) 
dt 



A^/3Ax^ a 
t" dx 

2t" 



dx 



dV{x,t) 



dx 



''nix, t) 



n{x,t) + 0{Ax^l2A) 



In the limit A — and r — > 0, with ~ 



2t° 



and 



ija = (2/3kq,)~^ we recover the fractional Fokkcr-Planck 
equation, Eq. ([?]), for subdiffusion in an external space- 
and-time-dependent force field 



F{x,t) = 



dV{x,t) 
dx 



(25) 



The fractional Fokker-Planck equation, Eq.(j3]), can 
also be derived from the subordinated stochastic 
Langevin equation motivated by physical arguments [loj 
to model subdiffusion in a space-and-time-dependent 
force field. This representation can be formulated as a 
system of stochastic equations. 



dYt 
dZt 



( F{Yt,Zt)r,-^ 







dt- 



{2K)^/^dBt 



dUt 



(26) 



where Bt is a one-dimensional Brownian motion and Ut 
is a a-stable Levy subordinator in [0,oo),0 < a < 1. 

I 



It is asumed that Bt and Ut are independent stochastic 
processes and the initial condition is Yq = Zq = 0. The 
stochastic process representing subdiffusion in a space- 
and-time-dependent force field is postulated to be given 
by = YiSt) where for t > 0, St is the random 
time the process Ut exceeds t. 

The stochastic differential equation, Ea. ((26l) . belongs 
to the general class of stochastic processes driven by Levy 
noise [20|, eq.(6.12)]. The infinitesimal generator for the 
process Eq. ([25|) is then given by 0, eq.(6.42)] 

Afiy,z) = ^^fiy,z) + n^Jiy,z) 



+ [f{y,z + z')-f{y,z)]—^z 
Jo r(i - a) 



'~^-"dz' 



(27) 

The Fokker-Planck evolution equation for the probability 
density qt{y, z) of the process {Yt, Zt) is given by 



d_ 

dt 



{y,z) = A'^qt{y,z) 



(28) 



where 



A^f{y,z) = K^f{y,z) - |- {p^f{y,z)^ 

-oD':fiy,z) 



(29) 



defines the operator adjoint to A. 

Now we relate the densities ptix) and qt{y,z) of the 
stochastic processes Xt and {Yt,Zt) respectively. We 
write u for a particular (random) path of the lat- 
ter process, and note that the coordinates at time t, 
{Yt[ijj), Zt{uj)) and Xt[ijj), are functions of w. For a fixed 
interval / we define the indicator function 



5i{x) 



1 if .T G / 

if .T ^ / 



(30) 



together with the auxiliary function 



H{t',uj,Az) = 



Si{Yt,{uj)) if Zt'^{io)<t<Zt'-{uj) + Az 
otherwise. 



r 



(31) 



With the above notation we can write 

pt{x)dx = {Si{Xt{uj))) (32) 

where the angle brackets represent an ensemble average 
over all paths w. Given that Xt can be interpreted as the 
Y coordinate of the last position of the Markov process 



(Y, Z) before it exits the set R x [0, t], we have 

5i{Xt{Lo)) = Y,H{t',io,AZt>{u:)), (33) 



t'>0 



since all summands equal zero except for t' ~ St, in which 
case Yt'{uj) = Xt{uj). The jumps Az = AZt'{uj) are 
a Poisson point process on (0, oo) whose characteristic 
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-l-a 



We can now 



measure has the density p^^j^^Az" 
combine Eas. ([32|) . (|33p and use the compensation formula 
in 



21|, XII (1.10)], to write 



Pt[x)dx = { / i/(t',L.,Az)^- -d^zdt' 

J I \Ja Jo r(l - a) 

(34) 

After integrating over Az and using Ea. ([5T|) the right 
hand side simphfies further to 



dt'5j{YAuj))S[o.t]{ZAuj)) ^f^'^^^^l" ) . (35) 
r(l - a) ' 



The ensemble average is evaluated using the probability 
density qtiy, z) so that 



Pt{x)dx ^ / dt' / dy / dz qt'{y,z)— 

/ Jo J I Jo r(l - a) 



and thus 



Pt{x) = / dt' o//-"gt,(a;,t), 
Jo 



(36) 



where ^I^ " denotes the Liouville fractional integral of 
order 1 — a acting on t. It also follows that 



oD'rytix) 



qt'{x,t)dt', 



—pt{x) = ^Dfqt.{x,t)dt'. 



(37) 



(38) 



We now solve Eqs.([28l), for QD^qt'{x,t) and subsi- 
tute this into Ea. ([38|) to obtain 

d f°° f Id 



rj dx 



-—q,,{x,t))dt' 



(39) 



where the jump probabilities are evaluated at the start 
of the waiting time prior to jumping. Note that this 
equation can be derived from the Montroll- Weiss CTRW 
formalism (see e.g., 0) if and only if the jumping prob- 
abilities are independent of time. Our inclusion of time 
dependence is an ad-hoc generalization for time depen- 
dent jumps. Using Laplace transform methods as above 
then leads to the discrete space evolution equation 

dn A d^^" 

+pi{x,+i,t)n,+i{t)} . (42) 

It follows from Eqs. (16) and (17) that this evolution 
equation is the evolution equation for the loss flux in 
the generalized mater equation approach. After taking 
the spatial continuum limit of Ea. (j42p with Boltzmann 
weighted jumping probabilities we have 



9?! 

'dt 
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dt 



d'^n(x,t) 1 d fdV(x,t) , , 
dx^ rjaadx \ dx 



(43) 

This provides an interpretation of the fractional Fokker- 
Planck equation 



dP 
'dt 



O^t 



92 19^, ; 

KaTT^ —F(X,t) 

dx'' rja dx 



P{x,t) (44) 



and finally, using Ea. ([37l 



as an equation for the loss flux in subdifFusion in a space- 
and-time-dcpendent force field 

It is straightforward jTsj to obtain numerical solutions 
of the discrete space evolution equations for subdifFusion 
in space-and-time-dependent force fields, Eqs. (|T7|) . ((42|) 
using an implicit time stepping method with the frac- 
tional derivatives approximated using the LI scheme (2^ . 
The CTRW models described here are also easy to sim- 
ulate using Monte Carlo methods 15 1. 
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g-^Pt{x) = K^^oA Pt{x) - { -^—oDt Pt{x) 



dx 



77 



(40) 

where we have used qoo{x, t) = and qo{x, t) = 5(q q){x, t) 
and assumed t > 0. Equation (UOl) recovers the fractional 
Fokker-Planck equation for space-and-time-dependent 
forces, Eq.Q. 

A different fractional Fokker-Planck equation can be 
obtained from the following generalization of the CTRW 
particle balance equation 



n,[t) = n,(0)$(i)+ / \pr{x,^i,t)n,^i{t) 



+pi{x^+i,t)n,+i{t)^i:{t - t) dt, (41) 
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